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Automatic reduction of four-loop bubbles * mit-CTP 3326 
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We give technical details about the computational strategy employed in a recently completed investigation of 
the four-loop QCD free energy. In particular, the reduction step from generic vacuum bubbles to master integrals 
is described from a practical viewpoint, for fully massive as well as QED-type integrals. 


1. Introduction 

Vacuum integrals, i.e. integrals without exter¬ 
nal momenta (often also called tadpoles or bub¬ 
bles) , constitute an important class of multi-loop 
Feynman integrals. While the perturbative ex¬ 
pansion of quantities like the free energy can be 
directly expressed in terms of vacuum integrals, 
they also serve as essential building blocks for 
many other computations, being the coefficient 
functions in asymptotic expansions of diagrams 
with external legs, and encoding the ultraviolet 
behavior of multi-scale integrals. 

A typical perturbative calculation proceeds in 
four conceptually independent steps. First, all 
relevant diagrams including their combinatoric 
factors are generated. For an algorithm that does 
this for vacuum integrals, see [1]. Second, the 
Feynman rules of the theory under consideration 
are inserted, and the color and Lorentz algebra is 
performed. Since in general individual loop inte¬ 
grals are divergent, a regularization scheme has to 
be adopted, the most practical one at present be¬ 
ing dimensional regularization (DR). Third, lin¬ 
ear relations between the regularized integrals are 
exploited, to systematically reduce all integrals 
occurring in the computation to a small set of 
so-called master integrals. In the framework of 
DR, the most important class of relations can be 
derived from integration-by-parts (IBP) identities 
[2]. Fourth, the master integrals have to be evalu¬ 
ated, either, in some fortunate cases, fully analyt¬ 
ically, or as an expansion in terms of the regular¬ 
ization parameter, in which case - and only here 
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- the number of dimensions d has to be specified. 
For results on the 4-loop level, see [3] (d = 4 — 2e) 
and [4] (d = 3 — 2e). 

At higher loop orders, it is inevitable to auto¬ 
mate the above setup to a large degree. There 
exist many approaches to implement automated 
perturbative calculations, and this is not the 
place to give a comprehensive review (see e.g. 
[5]). Instead, it is the third of the above steps 
that we wish to elaborate on in this contribution. 

A computer algebra system that is particularly 
well suited to cope with the demands of higher or¬ 
der perturbative calculations is FORM [6]. While 
by no means mandatory to use, we have adopted 
it to implement our algorithms, and hence we 
will indicate in a few places which specific FORM 
commands turned out to be extremely helpful. 

2. Notation and general considerations 

Consider the generic vacuum topologies of 
Fig. 1. In this intuitive graphical notation, ev¬ 
ery line represents a propagator (pf + mf)~ ai , 
with integer power at > 0, where the index i la¬ 
bels the different lines with momenta pi, which 
in turn can be expressed as a linear combination 
of the £ loop momenta kj. The vertices do not 
have any structure, except for assuring momen¬ 
tum conservation. Each diagram can carry a non¬ 
trivial numerator structure, which in the general 
case consists of powers of scalar products of the 
loop momenta. At £ loops, there are £{£ + l)/2 
different combinations ki • kj. 

Let us distinguish three different representa¬ 
tions of our integrals, which naturally appear at 
various levels of the reduction process: generic 
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Figure 1. The 1+1+3+10 generic vacuum topolo¬ 
gies up to four loops. The 0+1+2+6 factorized 
topologies are not shown here. 


integrals, their standard representations, and the 
master integrals. The goal of step three is then to 
formulate the algorithms which transform generic 
to standard to master integrals: 


ni<i<^ (fc* ■ kj ) b ' 3 
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Above, the label ’Filter’ symbolizes a collection 
of low-level routines, whose main action is to com¬ 
plete squares in the numerator and cancel against 
propagators such that only irreducible numera¬ 
tors remain. At this point, it is possible to repre¬ 
sent an Uloop vacuum integral by a list of t?(£ + l) 
non-negative numbers ({cq}, {bij, mj), the first 
half of them collecting the powers of propagators 
ai, while the second half contains either the power 
of an irreducible numerator (if the corresponding 
ai is zero) or the mass of the line. Furthermore, 
at this step equivalent topologies are re-labeled in 
a unique way by shifting the loop momenta, i.e. 
assigning a characteristic pattern of zeroes among 
the a,; to each topology of Fig. 1. 

In the remainder, we will specialize on two dif¬ 
ferent general classes of vacuum diagrams. First, 
we will consider all lines to have the same mass, 
mi = m. This class of integrals is useful when 
computing infrared-safe quantities like renormal¬ 
ization coefficients, in which case the infrared sec¬ 
tor of individual diagrams can be regulated by in¬ 
troducing masses into massless propagators. Sec¬ 
ond, we will allow for all ?n,; to be either zero or m, 


with the restriction that the number of massive 
lines at each vertex be even. This includes theo¬ 
ries like QED and gauge+Higgs models, whence 
we call this class ’QED-like’. 

The label ’Tables’ in Eq. (1) symbolizes a 
lookup in a database, which contains the neces¬ 
sary relations in a tabulated form. These tables 
are the main ingredient of the reduction step, and 
their organization and generation, which system¬ 
atically exploits IBP identities, will be described 
in more detail below. 

The intermediate step of applying the ’Filter’ 
algorithms not only serves the purpose of allowing 
for a fairly compact representation of the integral, 
but can also be used to keep the number of en¬ 
tries in the database, the memory requirements, 
and the CPU time needed for their derivation, in 
manageable bounds. To this end, we found it ad¬ 
vantageous to add further routines to the ’Filter’ 
package: 

• Early detection of zeroes: massless (sub-) 
tadpoles are zero in DR, as are integrals 
whose integrand does not depend on one of 
the loop momenta. 

• Symmetrization of the integrand: use the 
full symmetry group of the corresponding 
topology to order the list, and hence enable 
early cancellations in big expressions. 

• Decouple scalar products involving the loop 
momentum of a factorized one-loop tadpole: 

vanishes for odd n and is pro¬ 
portional to a totally symmetric combina¬ 
tion of metric tensors g^ lli2 ■ ■ ■ g^^}- 
The FORM function dd_ is perfectly suited 
for this symmetrization. This eliminates 
the need to derive relations for 8 (out of 
9, the 9th being the two-loop x two-loop 
case) of the factorized topologies, since af¬ 
ter decoupling the numerator, factorization 
into scalar vacuum integrals of the type of 
Eq. (1) is complete. 

• Reduce powers of factorized one-loop tad¬ 
poles to one: 

f d d k _ d — 2a f d d k 

J (k 2 + m 2 ) a + 1 2am 2 J (k 2 + m 2 ) a 

• Employ the ’triangle relation’ [2]: for in¬ 
tegrals involving massless lines, this helps 
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to reduce the number of different topolo¬ 
gies that have to be treated in the database 
considerably. 

Another potentially useful routine, which we have 
however not implemented, would be to use T- 
operators [7] in order to trade all numerator 
structure for higher dimensions of the integral 
measure, hence also immediately decoupling the 
factorized (two-loop x two-loop)-topology. 

One more practical note: To not miss cancella¬ 
tions, it is important to have a unique representa¬ 
tion for coefficients. Partial fractioning of terms 
like ~t ~— and —1— —A — helps here, ensuring the 
coefficients Cj(d) to be a sum of powers of simple 
poles and powers of d. 

In principle however, all these further relations 
are redundant since they would be automatically 
covered by the IBP identities. As mentioned 
above, their sole purpose is to optimize the deriva¬ 
tion of relations among the integrals, to be dis¬ 
cussed next. 

3. Reduction 

Integration by parts relies on the fact that an 
integral over a total derivative of any of the loop 
momenta vanishes in dimensional regularization. 
For the case of vacuum integrals, which we are 
interested in here, the IBP identities read 

0 = J {{ a i}> {bij, m i}) ) (2) 

where p,q € {fci,..., k/} cover all £ 2 different i- 
loop identities, and we have made use of the stan¬ 
dard representation introduced above. 

There are two possible general strategies imple¬ 
menting the IBP identities to find relations use¬ 
ful for reducing the integrals from their standard 
representation to master integrals. 

The first strategy is to derive general relations, 
valid for symbolic list-entries. These general sym¬ 
bolic relations can then be applied repeatedly to 
any integral of the specified class, no matter how 
large the powers are, to achieve the reduction. 
In practice however, it turns out that it is quite 
an art to shuffle IBP identities for integrals with 
symbolic indices such as to obtain useful reduc- 
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Figure 2. The 1+0+2+10 master integrals of 
QED type, up to four loops. Full lines carry a 
mass to, dotted lines are massless. All numera¬ 
tors are 1, all powers of propagators are 1. Note 
that there is no two-loop representative needed. 

tion relations. In absence of a generic algorith¬ 
mic formulation, it involves extensive handwork, 
and typically there are many special cases to be 
considered when pre-factors vanish at special pa¬ 
rameter values. At lower loop orders, there are 
complete solutions, see e.g. [9] for two-loop two- 
point functions with general masses, or [10] for 
three-loop vacuum integrals with one mass. 

The second strategy, nowadays constituting the 
mainstream of higher-loop computations, is a 
more brute-force approach, which however has 
the huge advantage of being perfectly suited to 
be completely automated. The main idea is to 
write down IBP identities for specific values of 
the indices. Introducing a lexicographic order¬ 
ing among the integrals [8], it is then possible 
to solve every single one of the IBP identities for 
the ’most difficult’ integral occurring. By starting 
from simple topologies (low number of lines), one 
systematically generates relations which express 
’difficult’ integrals in terms of ’simpler’ (in the 
sense of the ordering) ones. Solving an adequate 
set of fixed-index IBP relations, it is possible to 
express every integral of interest in terms of a few 
simple ones, ultimately the master integrals. 

The reason why the second strategy is suffi¬ 
cient for most computations is that in practice, 
one does not meet the most general integrals, but 
only a subset, typically characterized by an upper 
cutoff on the sum of indices. Indeed, dealing with 
a concrete model like QCD, knowledge of the ver¬ 
tex and propagator structure allows to constrain 
the set of possible indices ({a i}, {bij, mi}), hence 
rendering the search-space to be covered with IBP 
identities finite. 

Building up the relations proceeds as follows: 
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• Pick a list of indices ({cq}, {bjj. rrii}) that is 
’simple’, typically meaning a low number of 
loops, a low number of different lines, a low 
number of extra powers on the propagators, 
a low number of powers on the irreducible 
numerators. In FORM, these lists are most 
naturally represented by sparse tables. 

• Generate the first of the IBP identities. 

• Call ’Filter’ to transform the resulting sum 
of integrals to the standard representation. 

• Label the ’most difficult’ integral, according 
to the lexicographic ordering. A global op¬ 
eration like that became possible with the 
introduction of ’$-variables’ in FORM v3. 

• Invert its coefficient and multiply it into the 
equation. We do so only if we can factor¬ 
ize the coefficient into terms which are lin¬ 
ear in d, to preserve the generic structure of 
coefficients. While (at present) there is no 
factorization algorithm in FORM, we im¬ 
plemented one by ’guessing’ zeroes, utiliz¬ 
ing the fact that since the coefficients are 
generated by IBP, most of them have fac¬ 
tors ( nd ± a) where n is not bigger than the 
number of loops, and a is an integer of mod¬ 
erate size. To check that no relations are 
missed when factorization fails, it is useful 
to keep track of those cases and check in the 
end. 

• Bring the most difficult integral (having 
coefficient 1) to the left-hand-side, taking 
the generated equation as a definition. In 
FORM, this is done by the fill statement. 

• Take the next IBP identity, repeat the 
above steps. Increase list-indices. Repeat... 

• Write the relations found to disk in inter¬ 
vals. Large intervals ensure a high degree 
of re-substitution (of relations for integrals 
that are found later but that appeared on 
the right-hand-sides earlier), but are risky 
when the program execution crashes. 

Solving the IBP relations one by one like de¬ 
scribed above seems to be simpler than solving 
large systems of linear equations at once. In the 
end, it might be advantageous to re-substitute re¬ 
lations, which is possible by re-loading sets of re¬ 
lations into memory and re-writing them to disk. 
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Figure 3. The 0+0+0+3 fully massive master in¬ 
tegrals, in addition to those 1+1+3+10 of Fig. 1, 
taken at powers and numerators 1. A dot on a 
line means it carries an extra power. 


In the end, one has to check whether the 
set of generated identities is sufficiently large to 
achieve a reduction of all integrals occurring in 
the physics problem at hand. While a first edu¬ 
cated guess on the maximum powers needed can 
be obtained by scanning the terms to be calcu¬ 
lated after application of the ’Filter’ package, it 
might be necessary to enlarge the set of relations 
in further runs. To this end, the tablebase state¬ 
ment of FORM, implemented in version 3.1, al¬ 
lows for a good control over large amounts of data 
in the form of tables and table elements. 

The resulting master integrals are depicted in 
Fig. 2 for the ’QED-like’ case, and in Figs. 3,1 for 
the fully massive case. 

4. Master integrals 

Once the reduction algorithm ’stops’, are we 
guaranteed to arrive at the desired minimal set 
of master integrals? If we had followed the path 
of deriving generic reduction relations, valid for 
symbolic indices, the answer would be yes. For 
the implementation in terms of specific indices, 
one can however not be absolutely sure not to 
miss a relation which would only be detected 
when increasing the upper cutoff on indices of the 
integrand. For most practical purposes it might 
already be sufficient to work with an incomplete, 
but small, basis. 

In the case of gauge theories, it is also amusing 
to watch the gauge-parameter dependence as an 
indicator of how ’close’ one is to the minimal set, 
since in a full reduction gauge-parameter depen¬ 
dent terms cancel at an algebraic level, in d di¬ 
mensions, before evaluating the master integrals. 

The basis of master integrals is of course not 
unique, but depends on the actual choice of the 
lexicographic ordering. While we label an integral 
with unit numerator as ’simpler’ than one with in- 
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Figure 4. Relations for a basis conversion from the set of massive masters found in [3] to our notation. 


creased powers on the lines, one could as well have 
put priority on always reducing powers of propa¬ 
gators to one. The latter choice was adopted in 
a recent paper [3], where the three master inte¬ 
grals of Fig. 3 are replaced by ones of equivalent 
topology, but with irreducible numerators. Con¬ 
sequently, there must be linear relations between 
the two choices of basis, valid analytically in d di¬ 
mensions. From our tables, we simply read them 
off, see Fig. 4. 


5. Discussion 

We did not comment on the problem of so- 
called spurious poles here. Spurious poles are 
singular pre-factors, which can occur in the re¬ 
duction relations. They are difficult to avoid in 
general, if one is not willing to specify the dimen¬ 
sion yet in the reduction process. However, in our 
four-loop computation of the QCD free energy, we 
treated them after the reduction was performed 
successfully, by changing basis with the help of 
the tables. 

In principle, the package at hand can be used 
for other calculations requiring a four-loop reduc¬ 
tion of massive vacuum bubbles. One such ap¬ 
plication would be the re-evaluation of the QCD 
beta function. 


Acknowledgments: It is a pleasure to thank 
K. Kajantie, M. Laine and K. Rummukainen for 
an enjoyable collaboration on the matter pre¬ 
sented in the above. This work was supported 
in parts by the DOE, under Cooperative Agree¬ 
ment no. DF-FC02-94ER40818. 

REFERENCES 

1. K. Kajantie, M. Laine and Y. Schroder, Phys. 
Rev. D 65, 045008 (2002). 

2. K. G. Chetyrkin and F. V. Tkachov, Nucl. 
Phys. B 192, 159 (1981); F. V. Tkachov, 
Phys. Lett. B 100 , 65 (1981). 

3. S. Laporta, arXiv:hep-ph/0210336. 

4. A. Vuorinen, Master’s Thesis, Helsinki Uni¬ 
versity, 2001. 

5. R. Harlander and M. Steinhauser, Prog. Part. 
Nucl. Phys. 43, 167 (1999). 

6. J.A.M. Vermaseren, arXiv:math-ph/0010025; 
http://www.nikhef.nl/~form/. 

7. O. V. Tarasov, Phys. Rev. D 54, 6479 (1996). 

8. S. Laporta, Int. J. Mod. Phys. A 15, 5087 

( 2000 ). 

9. O. V. Tarasov, Nucl. Phys. B 502, 455 (1997). 

10. D.J. Broadlrurst, Z. Phys. C 54 (1992) 599; 
L. V. Avdeev, Comput. Phys. Commun. 98, 
15 (1996); P. A. Baikov and M. Steinhauser, 
Comput. Phys. Commun. 115, 161 (1998). 
















